Charged inclusion in nematic liquid crystals 
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We present a general theory of liquid crystals under inhomogeneous electric field in a Ginzburg- 
Landau scheme. The molecular orientation can be deformed by electric field when the dielectric 
tensor is orientation-dependent. We then investigate the influence of a charged particle on the 
orientation order in a nematic state. The director is aligned either along or perpendicular to the 
local electric field around the charge, depending on the sign of the dielectric anisotropy. The 
deformation becomes stronger with increasing the ratio Ze/R, where Ze is the charge and R is the 
radius of the particle. Numerical analysis shows the presence of defects around the particle for large 
Ze/R. They are nanometer-scale defects for microscopic ions. If the dielectric anisotropy is positive, 
a Saturn ring defect appears. If it is negative, a pair of point defects appear apart from the particle 
surface, each being connected to the surface by a disclination line segment. 

PACS numbers: 61.30.Dk, 61.30.Jf, 77.84.Nh, 61.30.Gd 
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I. INTRODUCTION 

Recently, a number of complex mesoscopic structures 
have been observed with addition of small particles into 
a liquid crystal matrix [l|jl3- I'^ nematics, inclusions can 
distort the orientation order over long distances, often 
inducing topological defects |l,|l,|a,B0l3l- As a result, 
anisotropic log-range interactions mediated by the dis- 
tortion are produced among the immersed particles. It 
also leads to the formation of various uncommon struc- 
tures or phases, such as chain aggregates [J, soft solids 
supported by a jammed cellular network of particles H, 
or a transparent phase including microemulsions |l(]llllj| . 

The short-range anchoring of the nematic molecules 
on the inclusions surface is usually taken as the origin 
of the long-range distortion Illlllll[lll[ll[ll. In 
the continuum approach in terms of the director n in a 
nematic state, the anchoring free energy is expressed as 
the surface integral. 



F„ 



Wa / dS{n ■ vY 



(1.1) 



where dS is the surface element, v is the normal unit 
vector at the surface, and Wa is a parameter represent- 
ing the strength of anchoring [Ij| . As another anchoring 
mechanism, electrically charged inclusions should induce 
alignment of the nematic order in their vicinity 15, 16] . 
This mechanism is relevant for ions which are naturally 
present or externally doped and for colloidal particles 
whose surfaces are highly charged. A large charged par- 
ticle can also be inserted into nematics. However, the 
effect of charges in liquid crystals remains poorly under- 
stood and has rarely been studied, despite its obvious 
fundamental and technological importance. It is of great 
interest how the charge anchoring mechanism works and 
how it is different from the usual short-range anchoring 
mechanism. 

It is also worth noting that the ion mobility in ne- 
matics is known to be anomalously low as compared to 
that in liquids with similar viscosity '17'|. Originally, de 



Gennes [I3 attributed the origin of this observation to 
a long-range deformation of the orientation order of the 
surrounding liquid crystal molecules. Analogously, small 
ions such as Na"*" or Cl~ in polar fluids are surrounded by 
a microscopic solvation shell composed of polar molecules 
aligned along the local electric field [I^ |20, |21| . 

The coupling of the electric field and the nematic ori- 
entation arises from the fact that the dielectric tensor 
Eij depends on the director n or the local orientation 
tensor Q.^ . The alignment along a homogeneous electric 
field is a well-known effect [13, while the alignment in 
an inhomogeneous electric field poses very complicated 
problems. We mention an experiment J2J|, in which an 
electric field was applied to nematics containing silicone 
oil particles to produce field-dependent defects. In this 
paper, we will study the orientation deformation around 
a single charged particle in a nematic state in the phe- 
nomenological Landau-de Gennes scheme using the ori- 



entation tensor Qij |22l l2.1 | . The optimal alignment min- 
imizes the sum of the Landau-de Gennes free energy and 
the electrostatic energy. Similar approaches have recently 
been used to calculate the solvation free energy of ions 
in near-critical polar fluids |2j . 

This paper is organized as follow. In Sec. II, we will 
present a general Ginzburg-Landau framework. By min- 
imizing the free energy functional, we will obtain equilib- 
rium equations satisfied by Qij , while the electric poten- 
tial obeys the Poisson equation with a dielectric tensor 
dependent on Q^. In Sec. Ill, we will estimate the free 
energy contributions around an isolated charged particle 
to find the range of the orientation deformation and the 
condition of strong deformation. Some discussion will 
be given on the anchoring effect around a charged col- 
loidal particle surrounded by counterions. We will also 
give some discussions for a charged colloidal particle sur- 
rounded by counterions. In Sec. IV, we will numerically 
solve the equilibrium equations. We shall find the for- 
mation of topological defects for large charge Ze and/or 
small radius R of the particle. 



II. GENERAL THEORETICAL BACKGROUND 

A. Free energy functional with charges 

The liquid crystal order is described in terms of the 
symmetric, traceless orientation tensor Qij ;17> ■22> ■23| . 
which may be defined as 



Qij{r) = {uiUj - Si.j/3), 



(2.1) 



where u denotes the unitary orientation vector of the 
nematic molecules. We introduce rotationally invariant 
quantities, 



J2 



E 



Q 



J3 



ijk 



^ij^jk^ki 



(2.2) 



We assume the Landau-de Gennes free energy Fi^ = Fo- 
Fg [l3| consisting of two parts. 



^n = dr 






ijk 



(2.3) 
(2.4) 



where V^ — d/dxk with xi — x, X2 — y, and x^ — z. 
Here ^ is a temperature-dependent constant and is posi- 
tive in the nematic phase, while B, C, and L are positive 
constants in the isotropic and nematic phases. The tran- 
sition is first-order for nonvanishing _B, and L ~ /csT/cr 
with a being the microscopic molecular length. The 
second part Fg is the gradient free energy representing 
an increase of the free energy due to inhomogeneity of 
Qij. We neglect another gradient term proportional to 

Next we consider the free energy contribution i^e aris- 
ing from the electrostatic interaction. It depends on the 
charge density p(r) and the polarization vec tor p(r ) of 
the liquid crystal molecules and is written as 



;or w(r) 

laiii 



F, 



dr 



2 E x'^PiP, 

ij 



1 



(2.5) 



where E = —Wcj) is the electric field with 0(r) being 
the electrostatic potential. The electric induction D = 

E + 47rP satisfies 



V • D = -V 



47rV • p = 47r/5, 



(2.6) 



where p{r) is the charge density. In our theory, it is 
crucial that the tensor x*-' in Eq.(2.5) depends on Qij (see 
Eqs.(2.11) and (2.13) below). However, we assume that 
it is independent of E^ neglecting the nonlinear dielectric 
effect I23 • We assume no externally applied electric field 
and set 



= 0, 
on the boundary walls of the container. 



(2.7) 



If infinitesimal space-dependent deviations 5p and 6p 
are superimposed on p and p, the incremental change of 
-Fe is given by 



<5F, - 



I dr ^ i^PiPjSx'^ +P]X''5pt 



-E-Sp^ 



where we have used the relation, 



-S{E'^) = ~\7{(t)SD) + iTT{(t)Sp - E ■ 5p). 



(2.8) 



(2.9) 



We minimize F^ with respect to p at fixed p and Qij. 
From SFc/Sp = we obtain 



E x'^Pi " ^» °^ P* == E X'l^i 



(2.10) 



Thus the inverse matrix of x''^' is the electric susceptibility 
tensor Xij related to the dielectric tensor Sij by 



- 47rxj, 



(2.11) 



The electric induction is of the usual form Di = ^ ■ SijEj 
and the Poisson equation (2.6) becomes 



V • U = - ^ V^EijVjCJ) = in p. 



(2.12) 



In this paper the dielectric tensor is assumed to be lin- 
early dependent on Qij as 



£.y(r) ^ so6ij +eiQi.j{r), 



(2.13) 



where eq is positive, but £1 is positive or negative depend- 
ing on the molecular structure [i3- Now, from Eq.(2.10), 
Fp is of the standard form, 



Fe^^j drY,e,,E,E,. 



(2.14) 



Y.U X'^X'^^XkE from Y.i Xik} 



,kj _ 



We obtain 6x^-' - 

6ij in the right-hand side of Eq.(2.8). Use of Sxij — 

eiSQtj/4:TT and Eq.(2.10) yields 



SF, 



dr 



J2 ^E,E,SQ, 



(2.15) 



The total free energy is the sum F = Fq + Fg + 
Fe. In equilibrium we impose the minimum condition 
5F/5Qij = X6ij, where A is the Lagrange multiplier 
ensuring the traceless condition J^i Qa — 0- Then 
Eqs.(2.3), (2.4) and (2.15) give 

(-A + CJ2)Q^J -b(Y1 Q^kQkj - \j2S.. 



£1 

Stt 



-LV'^Q^j - — EiEj - --E^^j p , 



(2.16) 



where we have ehminated A making all the terms trace- 
less. The field-induced change of Qij arises from the 
term on the right-hand side bilinear in E. It is nonva- 
nishing only for ei 7^ 0. In this work we are interested 
in the case of highly inhomogcncous E in the nematic 
phase. In the simplest case of weak, homogeneous E in 
the isotropic phase, Qij is simply given by the right-hand 
side of Eg. ( 2. 16) divided by \A\ to leading order in the 
field [ll|26|. 



B. Uniaxial and biaxial orientations 

We here diagonalize the orientation tensor Q. Let the 
maximum eigenvalue of Q be written as 25*1/3. Then the 
other two eigenvalues may be written as — S'i/3± 5*2 with 
Si > 1 5*2 1 . Use of the principal unit vectors, n, m, and 
£ = n X m, diagonalizing Q yields 



1 



Si{nn--I) + S2{mm - ££) 



(2.17) 



In this limit the dielectric tensor (2.13) is of the standard 
uniaxial form [131, 



e\\ninj + e±{Sij - niUj), 



(2.24) 



where e\i = £0 + 2£iS'/3 and e± — Eq ~ eiS'/3. From 
Eq.(2.15) an infinitesimal change of the director n -^ 
n + Sn induces a change in Fe given by 



5Fe^- [ dr — (E ■ n)E ■ Sn. 
J 47r 

Thus the equation for n is written as [isL Il6| 

eiS 



K{V^n)j 



in 



{E-n){E) 



0, 



(2.25) 



(2.26) 



where (• • •)^ denotes taking the perpendicular part of 
the vector (_L n). In the numerical analysis in this paper, 
however, we will not use the above equation in terms of n, 
because the defect-core structure can be better described 
in terms oiQij [2ll23l|. 



where / = {%} is the unit tensor. The orientation is 
uniaxial for 5*2 = but becomes biaxial for 5*2 9^ |17|. 
In terms of ^i and ^2 in Eq.(2.17), J2 and J3 in Eq.(2.2) 
are calculated as 



^2-3^? 



2^2', 



J.-'^sl 



Here, if we set 

5*1 = Scosx, 



5, = 3-1/25. 



Zo\02 • 



smx, 



we find 



J2 



J^ 



9 



^'cos(3x), 



(2.18) 



(2.19) 



(2.20) 



where the angle x is in the range — 7r/3 < X < 7r/3 from 
S'l > |52|. As is well-known, when Fq in Eq.(2.3) is 
minimized with i? > 0, the uniaxial orientation x = 
(or S2 = 0) is selected and S{= Si) is determined by 



2CS^ -BS-?,A = 0. 



Since 5 > we have S = B/AC +[{B/ACf + iA/2CY'^. 
Far below the transition temperature, the orientation 
is uniaxial and the amplitude S may be treated as a pos- 
itive constant outside the defect-core region. Mathemat- 
ically, this limit can be conveniently achieved if we take 
the limit of large A with B/A and C/A held fixed. The 
free energy i^o + Fg is then approximated by the Frank 
free energy, 

Fp^ Jdr^Kj2i^kn^)^ (2.22) 

ki 



with the single Frank coefficient, 
K = 2S^L. 



(2.23) 



C. Axisymmetric orientation 

As an application of our general framework, we focus 
on the basic problem of a single charged particle im- 
mersed in a nematic liquid crystal. The electric field 
around the particle induces a deformation of the orienta- 
tion order parameter. Far from the particle the orienta- 
tion is assumed to be uniaxially along the z-axis. Before 
treating this case we here present the equations for Qij 
for general axisymmetric orientation. 

It is convenient to use the cylindrical coordinates 
(p, (fi, z), where p = (x^ -I- y^)^'"^ and tan(/3 = y/x. Using 
the unit vectors Sp — {x/ p,y/ p,Q), e^p — (— y/p, x/p, 0), 
and Ez = (0,0,1), we express the traceless orientation 
tensor as 

Q = Qiepep + Q2ezez + Q3{epez + e^Bp) 

~ (Qi+Q2)e^e^, (2.27) 



(2.21) where Qi {i — 1, 2, 3) depend on p and z. The eigenvalues 



of Q are given by A± = (Qi + Q2)/2 ± [Qj + (Qi 



Qi ^ Q2- Here we limit ourselves 



74] 1/2 and A, 
to the case A+ > X^ or 

3(Qi + Q2) + [4Qi + (Qi - Q2)']'/' > 0, (2.28) 
under which Si and S2 in Eq.(2.17) are written as 

3, 



^1 
^2 = 



4Ql + (Qi-Q2)']^/'- 



(Q1 + Q2), (2.29) 



ml + ((Qi - Q2)']'/' - 7(Qi + Q2). (2.30) 



With these definitions we shall find 5*1 > 5*2 > every- 
where in our numerical calculations, where the equality 
^i = S2 holds in an exceptional case (see Eq.(4.10)). The 



director n is perpendicular to e„ and the two components 



Bp ■ n and Uz — e^ ■ n satisfy 



[(Ql - Q2r + ^QiV^' + Qi-Q2}t^. (2.31) 

nz I ) ZQs 

In terms of Si, 5*2, rip, and Uz, we may express Qi as 

Ql = (^1 + 52)n2 - i^i - 52, (2.32) 



Q2 = 


= (5i + 52)n2 - -5i - 


-S2, 


(2.33) 


Qa = 


= iSi + S2)npnz. 




(2.34) 



In particular, Qi + Q2 = 5i/3 — 52- In our numerical 
analysis to follow, the condition (2.28) will be satisfied. 
However, if the reverse relation of Eq.(2.28) holds, we 
have 25*1/3 = — Ql — Q2 > and n — e^. 
From Eq.(2.19) J2 and J3 are written as 



J2 = 2{Ql + Ql + Ql + QiQ2), 
J3 = 3(Qi + Q2)(Q^-QiQ2). 



(2.35) 



The Fq can then be expressed in terms of Qj. The gra- 
dient free energy Fg in Eq.(2.4) is of the form, 



Fg = L f dr 


(VQi)^ + (VQ2)2 + (VQ 


+ WQi-\ 


,^ , (2Qi + Q2)2 , Ql] 

'Q^+ p2 + p2\ 



^ , (2.36) 
P- 

where (VQj)^ = {dQ^/dzf + {dQ.Jdpf and the last two 
terms in the integrand (oc p^^) arise from the relations 
dep/dip = e^p and de^jd^p = —Bp. 

Also the electric potential = 0(/9, z) is a function of 
p and 2; and the electric field is expressed as 



E^Eiep + E2ez 



(2.37) 



where Ei — — 90/9p and E2 = —d(j>/dz. The electric 
induction D — *e ■ E is expressed as 



D = [(eo + eiQi)Ei+eiQ3E2]ep 
+ [{eo + eiQ2)E2+siQ3Ei]ez. 



(2.38) 



The Poisson equation (2.12) becomes 


1 d 

-po-p' 


(eo..iQi)|..iQ3f 


d ' 

^0-z_ 


(eo+eiQ2)^^+eiQ3^/l 
oz dp 



-Airp. (2.39) 



We now set up the equilibrium equations for Qi requir- 
ing SF/6Qi = 0. From Eq.(2.15) we derive the relation. 



6F,^~^ I dr{El6Qi+El5Q2 + 2EiE25Qz). (2, 



40) 



Some calculations yield 
AL 



D2 - iV" 



2 ,{2Qi + Q2)-D3^^El (2.41) 



87r 



2i 

P 
L 



I?2 - iV^ + ^ ) (2Q2 + Ql) - A3 = i^E'2, (2.42) 






D2-LV' + ^ \Q^ = -LEiE; 



£1 



(2.43) 



The left-hand sides are 5Fi,/5Qi with Fi, = Fq + Eg and 
the right-hand sides are —5Fe/SQi. We introduce 

D2 = -A + CJ2-B{Qi + Q2), 

D3 = B(2Q2 + 2Q2' + 5QiQ2-Q3'). (2.44) 

Here we may define the amplitude 5 = (3J2/2)-'^/^ and 
the angle x as in Eq.(2.19). Then some calculations give 

3V3[4Ql + (Ql - Q2)']'/'i^3 = 255^ sin(3x) (2.45) 

Thus D3 oc {dJs/dxjs, so D3 = leads to the uniaxial 
orientation 52 = 0. Further requirement of D2 = is to 
impose Eq.(2.21). 



III. ESTIMATIONS OF THE FREE ENERGY 
A. An isolated charged particle 

We estimate the free energy contributions around an 
isolated charged particle with charge Ze and radius R 
deeply in the nematic state. 

If ei = 0, no orientation disturbance is induced and 
the electric potential is given by Ze/sor, where r is the 
distance from the particle center. In this case, since 
E^ = (Ze/eo)^?'^'*, the electrostatic free energy Fe is 
dependent on the lower cut-off radius R as 



Ft, 



dr^E' 



r>R 



2eoR 



(3.1) 



This is the first theoretical expression for the solvation 
free energy of ions in a polar fluid, where the lower bound 
R is called the Born radius piil l2(l l2l| . For not large £1 
we may expand E,, as 



F, = ^R 



drY.^E.,E,Q,,+0{el), (3.2) 



which follows from Eq.(2.15) at fixed charge density. For 
the case |ei| <£o we use the above expansion with E = 
— {Ze/sor^)r to estimate Fg. Further assuming that the 
orientation is nearly uniaxial, we obtain 



F,-F^^-e,S—^ dr-{n-rf 

°~£o JR<r<i ' 



(3.3) 



For £1 > (for £1 < 0), n tends to be parallel (perpen- 
dicular) to f — r~^r near the charged particle R < r < £, 



while (n-r)'^ should be replaced by the angle average 1/3 
far from the particle r > £. We assume that the orienta- 
tion disturbance is strong in the region R < r < £. On 
the other hand, the Frank free e nerg y Fp in Eq.(2.22) 
is roughly of order 2'KK{i - R) [23. For ei > the 
change of the total free energy i^ due to the orientation 
deformation is estimated as [isl [lq | 



2^2 



AF 



2£i SZ^e 
3 2el 



]-h+2.Kii- 



R). (3.4) 



For £1 < the factor 2ei/3 of the first term should be 
replaced by |£i|/3. However, the numerical factors of the 
two terms in Eq.(3.4) are rough estimates and should not 
be taken too seriously. 

Minimization of the first two terms on the right-hand 
side of Eq.(3.4) with respect to £ gives 



= Ze{\ei\S/6TrelK)^/'^. 



(3.5) 



Strong orientation deformation occurs for £ > R, which 
may be called the strong solvation condition for an iso- 
lated charged particle in liquid crystals. To characterize 
the strength of the charge we hereafter use the length, 



Ae == Z{£BkBT/8TTK) 



1/2 



(3.6) 



where £b = e^ /e^kT is the Bjerrum length usually of 
order of 10 nm for liquid crystals. In terms of Ae the 
strong solvation condition £ > Ris written as 



(|£l|/£0)l/2A,/i?>l 



(3.7) 



for 5 ~ 1. The left-hand side of Eq.(3.7) represents the 
dimensionless strength of the charge-induced deforma- 
tion iljj. In the reverse case (|£i|/£o)^/^Ae/i? < 1, the 
orientation deformation is weak, which may be called the 
weak solvation condition. 

From Eq.(3.4) the minimum (equilibrium) value of Ai^ 
is written as 



AF = 2tiK{2£ - ^V^ - R). 



(3.8) 



in terms of £ in Eq.(3.5). Here the right-hand side is 
zero ioT £ ^ R and decreases for larger £. Of course, 
even if Eq.(3.7) does not hold, weak deformations are 
induced to make AF < as long as £1 7^ 0. Notice that 
we are neglecting such weak deformations in the present 
estimations. See Fig. 14 below for numerical data of AF. 



Z/R, for example, when the ionizable points on the sur- 
face is proportional to the surface area AttR^. However, 
the problem can be very complex, because the small 
countcrions can induce large deformations of the orien- 
tation order around themselves. 

For simplicity, let us assume that the counterions do 
not satisfy Eq.(3.7) and only weakly disturb the orienta- 
tion order and that the screening length Ag of the colloid 
charge is shorter than R. Then the distribution of the 
counterions is close to that near a planar charged surface 
and As is given by the Gouy-Chapman length 28], 



A, = 2kBT/eE, = l/27r^Bff. 



(3.9) 



where ct^ = Z/AttR'^ is the surface charge density in units 
of e and Eg = AncTse/eo is the electric field at the surface. 
To ensure the inequality As < i? we require 



2TT£BRcr, > 1. 



(3.10) 



Then Fe — Fb in Eq.(3.3) may be written in the same 
form as Fa in Eq.(l.l), 



K = -\WejdS{n-rf. (3.11) 



Using Eq.(3.9) we obtain H 

We 



= SsiZ'^e^K/ATTelR^ 

= 2kBTSasei/eo, 



(3.12) 



which is independent of R if cts is a constant or Z oc R^. 
The anchoring due to the surface charge becomes strong 
for \We\R/K > 1 with K being the Frank constant in 
Eq.(2.23). This is analogous to the well-known strong 
anchoring condition |Wa|i?/iir > 1 for the neutral case 



IV. NUMERICAL CALCULATIONS 

For an isolated spherical particle without the counteri- 
ons, we numerically solved Eqs.(2.41)-(2.43) satisfied by 
the three components Qi and the Poisson equation (2.39) 
for the electric potential </>, without the microscopic an- 
choring interaction in Eq.(l.l). The particle radius is 
assumed to be considerably larger than the defect core 
radius. 



B. A charged colloidal particle 



Method 



Although it is not clear whether or not ionization can 
be effectively induced on colloid surfaces in liquid crys- 
tals, we here assume the presence of charged colloidal 
paticles in nematics. In such situations, the distortion of 
n due to the surface charge can be more important than 
that due to the anchoring interaction given in Eq.(l.l). 
The strong solvation condition (3.7) is satisfied for large 



To calculate the equilibrium Qi and (f>, we integrated 
the time-evolution equations. 



^Q^ 

dt 
9t 



5F_ 
-V ■ D + inp. 



(4.1) 
(4.2) 



See Eqs.(2.41)-(2.43) and the sentence below them for 
5F/5Q, and the left-hand side of Eq. (2.39) for V-D. The 
steady solutions reached at long times are the equilibrium 
or metastable solutions. In the following figures we will 
show the steady solutions only. 

We used a discretized 200 x 200 cell in the (/?, z) semi- 
plane (0 < p < 200Aa; and < z < 200Aa;) assuming 
the symmetry around the z axis and with respect to the 
xy plane (see the comment (i) in the last section) . In Fl 
we set B/A = 1 and C /A = 3. We took the mesh size of 
the grid at 



Ax = {L/Ay^'/2 ^ ^/2 



(4.3) 



We will measure space in units of Aa; in the following 
figures. The length ^ — (L/A)^/^ is the shortest one in 
our theory and is of the order of the defect core size. The 
particle radius R was fixed as a relatively large value, 



i? = IOC = 20Aa;. 



(4.4) 



All the figures to follow will be given in the region —100 < 
X < 100 and —100 < z < 100 with the x and z axes being 
horizontal and vertical, respectively. 

On all the cell boundaries, we assumed the constant- 
potential condition cj) ^ in Eq.(2.7) and the uniaxial 
orientation along the z axis. 



1 



-Sh, Q2 — TT'S'b, Qz — 0. 



(4.5) 



Here Sh = 0.79 is the bulk amplitude obtained as the 
solution of Eq. (2.21). We calculated Qi only outside the 
sphere r > i?. If we superimpose SQtj on Qij with SQij ~ 
on the system boundaries, the incremental change of F 
is written as 



SF ^-L 



/ da^{u-V Qij) SQij, 



(4.6) 



47r/5 for r < i? and Eq.(2.39) for r > R. As the steady 
solutions of Eq.(4.2), we confirm that (p and r • D change 
continuously across the interface even on the scale of the 
mesh size Ax, which are required in electrostatics. 

Since R/^ is fixed at 10, the remaining relevant con- 
trol parameters are the ratios Xe/R and ei/eo- We thus 
varied Ag at Ei/eq — ±1. 




FIG. 1: Q2 and Q3 for (a) Xe/R = 10 and (b) 5 at ei/so = 1 
in the range -100 < a; < 100 and -100 < z < 100 in units 
of ^/2 (see Eq.(3.2)). The z axis is vertical and the x axis is 
horizontal. There is a Saturn-ring defect in (a) and no defect 
in (b). 



B. Orientation for ei > 



where bulk contribution vanishes, da is the surface ele- 
ment, and V is the outward normal unit vector on the 
particle surface (equal to r~^r here). To ensure 5F = 
we thus imposed 



V ■ VQy = 0, 



(4.7) 



at r = i?. This becomes r-VQi — pdQi/dp + zdQi/dz = 
at r = _R in the axisymmetric case in Eq.(2.27). We 
assume no surface free energy or no supplementary an- 
choring on the particle surface. That is, we set Wa = in 
Eq.(l.l). On the other hand, the electric potential was 
calculated in the whole region in the cell for the com- 
putational convenience. That is, we assumed isotropic 
polarizability [D = e^E) inside the sphere r < R and 
used the smooth charge-density profile, 

p(r) = Ipo - ^/5o tanh[(r - R)/^] , (4.8) 

in the whole region. The constant po is determined from 
/ drp{r) = Ze. This means that (j) obeys — eoV^0 = 



We first focus on the case ei/eo = 1, where the parallel 
alignment n || ^ is favored near the particle. This corre- 
sponds to the homeotropic alignment (perpendicular to 
the surface) in the neutral case. In Fig.l, we display Q2 
and Q3 in Eq.(2.27) in the xz plane for (a) Xe/R = 5 and 
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FIG. 2: Nematic director n in the xz plane for (a) X^/R = 10 
and (b) 5 at ei/eo = 1. A Saturn ring is present in (a) and is 
nonexistent in (b). 




FIG. 3; Landau-de Gennes free energy density fg + fo — fox 
for Xe/R ~ 10 and ei/eo = 1 in the xz plane, where /ooo 1 
the background value of /o. It is peaked at the Saturn ring. 
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FIG. 4: Scalar order parameter Si and biaxiality parameter 
S2 for \e/R ~ 10 and ei/eo = 1 in the xz plane. 



can see a -1/2 defect on both sides of the particle on the 
X axis for Xe/R — 10 in (a), but there is no defect for 
Ae/i? = 5in (b). 

In Fig. 3, we plot the deviation of the Landau-de 
Gennes free energy density fg + fo — /ooo at Xe/R = 10, 
where Fq = J drf^ and Fg = J drfg are given in 
Eqs.(2.3) and (2.4). The /ooo is the value of /o far from 
the particle in the uniaxial state with Si — Sb and S2 — 
(see Eq.(4.1)). This figure clearly demonstrates the ex- 
istence of a Saturn ring at the sharp peaks. Next, in 
Fig.4, we display ^i in Eq.(2.29) and S2 in Eq.(2.30), 
respectively, for Xe/R = 10. At the defect positions. Si 
becomes small and S2 exhibits a peak, while they tend 
to their bulk values. Si — 0.79 and S2 = 0, far from the 
defect. In the defect region, the nematic is locally melt 
and biaxial. The latter biaxial behavior was also found in 
the previous molecular dynamics simulation for a neutral 
particle under the homeotropic anchoring [j, Q . 

We performed simulations for various other parameters 
(not shown here). With increasing Xe/R, a Saturn ring 
appears suddenly with a nonvanishing radius at a certain 
transition value of Xe/R (see Fig. 9 below). Essentially 
the same behavior can be observed with increasing Si/sq 
at fixed Xe/R- At small Xe/R, the couphng between the 
electric field and the nematic order is not strong enough 
to induce the radial orientation on the equatorial plane. 
We therefore observe the quadrupolar symmetry without 
defect formation for Xe/R = 5 in Figs.l and 2, where Si, 
S2 and /o + fg exhibit no peaks around the particle. 




a) Q, 



b) Q, 



(b) 10. These quantities are related to n as in Eqs.(2.32) 
and (2.33). Figure 2 displays the corresponding config- 
uration of n = npBp + HzBz in the xz plane deduced 
from Eq.(2.31). For large Xe/R, the nematic exhibits a 
radial orientation at the surface of the particle and a line 
of —1/2 defect appears surrounding the particle in the 
plane z = 0, because the orientation is along the z axis 
far from the particle. Such a circular defect line is called 
"Saturn ring" in the literature [j,|j,l3. In Fig. 2, we 



FIG. 5: Q2 and Q3 for (a) Xe/R = 10 and (b) 4.2 at ei/eo 
— 1 in the xz plane. 



C. Orientation for ei < 

We turn to the second case ei < at ei/eo = —1, 
where the perpendicular alignment n _L £^ is preferred 



•. W I / ^ J' 






a) 



b) 



FIG. 6; Nematic director n for (a) Xe/R = 10 and (b) 4.2 at 
£i/£o = — 1 in the xz plane. A pair of point defects can be 
seen in (a), while there is no defect in (b). 




FIG. 7: Landau-de Gennes free energy density ,fo + fg — fooo 
for Xe/R = 10 and ei/eo ~ —1 in the xz plane, where fooo is 
the background value of /o. It is peaked along the disclination 
line segments connecting the poles and the point defects. 



close to the particle. We show Q2 and Q3 in Fig. 5 and n 
in Fig. 6 for (a) Xe/R = 10 and (b) 4.2 in the xz plane. 
Remarkably, in (a) the director is everwhere parallel to 
the surface and a pair of point defects are created on 
the z axis, while in (b) the director is not parallel to 
the surface in the neighborhood of the poles. With these 
defects in (a) , the Landau-de Gennes free energy density 
/o + fg exhibits twin maxima at the point-defects as in 
Fig. 7, while 5*1 and S2 behave as in Fig. 8. 

The defect structure in (a) has never been reported in 
the literature. Each point defect is detached from the 
surface and is connected to one of the poles by a +1 
disclination line segment with length £d. Here id — 37 
in (a). In fact, on both sides, a ridge structure with its 
top at the point-defect can be seen in /o + /g, S'l, and ^2 
in Figs. 8 and 9. The singular line segments created are 
specified \iy R < \z\ < R + id and p = 0, on which our 
numerical analysis yields Q3 = and Q\ — — Q2/2 > 0, 
so that 



(3 = (3i(/ -Se^e^). 



(4.9) 




Here / is the unit tensor. This form is natural from the 
rotational invariance around the z axis and n, = 0. On 



100 



FIG. 8: Scalar order parameter S\ and biaxiality parameter 
& for Xe/R — 10 and ei/eo = —1 in the xz plane. They 
deviate from the background values around the disclination 
line segments. 



these line segments, the equality, 

^1 = ^2 = 3Q1/2, (4.10) 

can be found from Eqs.(2.29) and (2.30). To see the 
defect structure in more dtail, we plot Si and iS'2 as func- 
tions of z for p — (along the z axis) and p = 3 in Fig. 9 
and S'l, 5*2, and nz as functions of p for 2: = 30 in Fig. 10. 
We can see how S'l becomes small and equal to S2 along 
the z axis, while Si > S2 not on the z axis. We also 
recognize that Uz changes slowly on the scale of R, while 
Si and S2 change on the scale of the core radius. A high 
degree of biaxiality is present around the line segments. 
In particular. Si = S2 on the singular line segments and 
Si > S2 outside them. 

In the neutral case with short-range anchoring, a pair 
of point defects, called boojums Jd^,2M' ^^^ attached to 
the surface at the poles for large negative Wa- 



D. Defect - no defect transition 

To gain quantitative information, we computed the dis- 
tance of the defect core from the particle center id as a 
function of Ae/i? for various values of Ei/eq. As shown in 
Fig. 11, a Saturn ring appears discontinuously for positive 
El . The solid curves were obtained when Ae was increased 
from small values, while the dotted curves were obtained 
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FIG. 9: 5*1 (solid line) and 5*2 (dotted line) along the z axis 
with Xe/R = 10 and £i/eo = —1 for (a) p = and (b) p = 3, 
where Si = S2 on the Une segment R < z < R + id at p = 0. 
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FIG. 11: Normalized defect radius (.djR vs Ae/-R for vari- 
ous ei/eo- For positive £1 the defect formation and collapse 
are discontinuous, where the solid and dotted curves are the 
results of increasing and decreasing Ae, respectively. For neg- 
ative £1 the transition is continuous with a critical point. 





FIG. 10; 5*1 (dotted line), ^2 (dotted line), and n^ (solid line) 
at z = 30 along the p axis for \e/ R = 10 and £i/£o = —1. 



FIG. 12: Normalized changes of the total free energy Eq.(4.7) 
(solid lines) and the Landau-de Gennes free energy Eq.(4.6) 
(broken lines) vs \e/ R for three £i/£o. For £1 > the discon- 
tinuities correspond to the defect formation in Fig.ll. 



when Ae was decreased from large values. We can see that 
the value of \e/R at the formation is larger than that at 
the disappearance. Similar hysteresis behavior was also 
found in a two-dimensional simulation with positive Wa 
by YamamotoQ. For negative ei, a pair of defects ap- 
pear from the surface continuously and is detached for 
Xe/R larger than the critical value (~ 4.9 in Fig.ll). 
For both positive and negative £1, the transition value 
of Xe/R decreases with increasing |ei|/eo (not shown for 
negative £1). These transitions are consistent with the 
criterion (3.7) since the defect formation is possible only 
in the strong solvation condition. Furthermore, once the 
defect is created, the defect size Id grows linearly with 
Xe/R, with the slope increasing with increasing |£i|/£o 
(which is the case also for negative £1). This trend is 
also consistent with Eq.(3.7), provided that the distance 
of the orientation deformation £ there and the defect size 
£d here are assumed to be of the same order. 



Next, in Fig. 12, we plot the increase of the Landau-de 
Gennes free energy, 



AFl= / dr[fo + fg - fo^], 

Jr>R 

and that of the total free energy, 



AF = F- F„„ - AFr 



AFe 



(4.11) 



(4.12) 



as a function of Xe/R for three values of £i/£o, where AFe 
is the increase of the electrostatic part. These quantities 
were calculated outside the particle r > R using the so- 
lutions of Eq.(2.39) and Eqs.(2.41)-(2.43). The Fu„ is 
the value of F in the unperturbed state in which the 
orientation is uniaxially along the z axis and the electric 
potential ((> is calculated using the homogeneous Sij in the 
uniaxial state. Thus these increases vanish as Xe/R -^ 0. 
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We recognize that AFl increases but AF decreases with 
increasing Xe/R- These aspects are consistent with the 
estimations in Eqs.(3.4) and (3.8), though they are rough 
approximations. It is remarkable that AF^ is negative 
and its absolute value is larger than AFj^. In addition, 
at the Saturn-ring formation, AF/LR changes by —1.1 
and —0.9 for ei/eo = 1 and 2, respectively, while the 
corresponding changes of AFi^/LR are 8.2 and 7.1, re- 
spectively. Notice also that AFl and AF are continuous 
even at the defect formation for negative ei. 



V. SUMMARY AND DISCUSSIONS 

We summarize our main results. In Sec. 2, we have 
derived the equilibrium equations for the tensor order 
parameter Qij as in Eq.(2.16), supplemented with the 
Poisson equation (2.12). The interaction between the 
orientation and the electric field arises from the dielectric 
anisotropy ei in Eq.(2.13). In the axisymmetric case 
Qij is expressed in terms of the three components Qi 
{i — 1,2,3) as in Eq.(2.27). They are determined by 
Eqs.(2.41)-(2.43) with the aid of the Poisson equation 
Eq.(2.39). In Sec. 3, we have estimated the range € of 
the strong orientation deformation around an isolated 
charged particle as in Eq.(3.5) and obtained the criterion 
of the strong deformation as in Eq.(3.7). There, we 
can see analogy between the present problem of liquid 
crystals and that of the ion solvation in polar fluids 2V^ . 
For a highly charged colloidal particle, we have found the 
effective anchoring parameter Eq.(3.12) under the con- 
dition Eq.(3.9). In Sec. 4, we have numerically examined 
the orientation deformation around an isolated charged 
particle. For positive £i and for large Xe/R oc Z/R, we 
have obtained a Saturn ring around the particle, where 
Ae is the characteristic length in Eq.(3.6). As shown in 
Fig.ll, the defect radius exhibits a discontinuous change 
as a function of Xe/R for positive £i. This means that a 
Saturn ring appears or disappears suddenly with radius 
larger than R as Xe/R is increased or decreased. For 
negative ei and large Xe/R, we have found appearance 
of point defects (boojums) on both sides of the particle 
along the z axis as a continous transition. They are 
detached from the surface, while they are on the surface 
for a neutral particle. As shown in Fig. 12, the decrease of 
the electrostatic free energy overcomes the increase of the 
Landau-de Gennes free energy, resulting in the lowering 
of the total free energy, with the orientation deformation. 

We make further comments on the limitations of our 
work and future problems. 

(i) In our simulations we have sought only the axisym- 
metric defrmations symmetric with respect to the xy 
plane. However, as in the neutral case |l,ll>IEIy|i there 
can be symmetry-breaking defects such as a combination 



of a radial and a hyperbolic hedgedog. 
(ii) The orientation deformation occurs much longer than 
the particle radius. For microscopic ions the deforma- 
tion extends over a nanometer scale. From Fig. 12 we can 
see that the lowering of the free energy much exceeds 
LR ^ ksTR/a, where a is the size of a liquid crystal 
molecule. Thus, even for i? ~ a, we find |AF| ^ ksT 
and the deformation can be stable against thermal fluc- 
tuations. As mentioned in Seel, the long-range deforma- 
tion should be the origin of anomalously low ion mobility 
in nematics [i3ll3- However, the existence of defects on 
a nanometer scale is not well established because of the 
limitation of our coarse-grained approach. 
(iii) We have examined the charge effect only in nematics. 
Slightly above the (weakly first-order) isotropic-nematic 
phase transition, the deformation around charged parti- 
cles can be long-ranged, extending over the correlation 
length '3Q3, as in near-critical polar fluids f21|. Fur- 
thermore, metastability or hysteresis behavior at the 
isotropic-nematic transition should almost disappear in 
the presence of a small amount of ions, owing to ion- 
induced nucleation, as in polar fluids J3lj . 
(iv) We have presented numerical analysis only for a sin- 
gle particle, but long-range correlations among charged 
inclusions can produce a number of unexplored efi^ects. In 
particular, doping of a small amount of ions into liquid 
crystals should provide highly correlated ionic systems, 
where transparent nematic states would be realized [lO| . 
Ion solubility and ion distribution across an isotropic- 
nematic interface should constitute new problems, which 
have been studied for polar fluids |3^ . 
(v) We should investigate dynamical properties of 
charged particles in liquid crystals such as the ion mo- 
bility or convection [13, ISfl ■ They should be much influ- 
enced by the orientation deformation jSJ] . 
(vi) Beyond the particular interest in the field of liquid 
crystals, the general approach developed in this paper 
and preceding ones [l^. Ha yB HB H^ using inhomoge- 
neous dielectric constant or tensor could provide a prac- 
tical and coherent method to study various polarization 
effects in simple and complex fiuids. 
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